packages<-c('ggplot2','dplyr','mada')
lapply(packages, library, character.only=T)
## [[1]]
## [1] "ggplot2"   "stats"     "graphics"  "grDevices" "utils"     "datasets" 
## [7] "methods"   "base"     
## 
## [[2]]
## [1] "dplyr"     "ggplot2"   "stats"     "graphics"  "grDevices" "utils"    
## [7] "datasets"  "methods"   "base"     
## 
## [[3]]
##  [1] "mada"      "mvmeta"    "ellipse"   "mvtnorm"   "dplyr"    
##  [6] "ggplot2"   "stats"     "graphics"  "grDevices" "utils"    
## [11] "datasets"  "methods"   "base"
"%&%" = function(a,b) paste(a,b,sep="")
h2.dir <- '/Volumes/im-lab/nas40t2/Data/Annotations/heritability/'
gcta <- read.table('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/expArch_DGN-WB_imputedGTs/DGN-WB.h2.all.models_FHSfdr0.05.Chr1-22_globaleQTLOtherChr_reml-no-constrain.2015-12-13.txt',header=TRUE) %>% dplyr::select(gene,local.h2)
dim(gcta)
## [1] 12719     2
##get price h2 (custom-made human array containing 23,720 unique oligonucleotide probes, http://www.nature.com/nature/journal/v452/n7186/full/nature06758.html)
price <- read.table(h2.dir %&% "Alkes/h2all.txt",header=TRUE) %>% dplyr::rename(gene=gname)
dim(price)
## [1] 15078     7
wright <- read.table(h2.dir %&% "Fred/h2PB.txt",header=TRUE) %>% dplyr::rename(gene=gname)
dim(wright)
## [1] 17792     2
all<-full_join(gcta,price,by='gene')
all<-full_join(all,wright,by='gene')
#add GTEx whole blood
gtex <- read.table('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/GTEx.TW.WholeBlood.h2.all.models_FHSfdr0.05.Chr1-22_globaleQTLOtherChr_reml-no-constrain.2015-12-14.txt',header=T) %>% mutate(gtexWB=local.h2) %>% dplyr::select(gene,gtexWB)
all<-full_join(all,gtex,by='gene')
dim(all)
## [1] 22567    10
nona <- all[complete.cases(all),]
dim(nona)
## [1] 5518   10
#dgn v. gtex whole blood
ggplot(nona,aes(x=local.h2,y=gtexWB))+geom_point(alpha=1/4)+geom_smooth()+theme_bw()+coord_cartesian(ylim=c(-0.15,1.02))

res<-cor.test(nona$local.h2,nona$gtexWB,method='s')
CIrho(res$estimate,length(nona$local.h2))
##            rho     2.5 %    97.5 %
## [1,] 0.2713322 0.2467124 0.2956019
#dgn v. price blood
ggplot(nona,aes(x=local.h2,y=h2bloodcis))+geom_point(alpha=1/4)+geom_smooth()+theme_bw()

res<-cor.test(nona$local.h2,nona$h2bloodcis,method='s')
CIrho(res$estimate,length(nona$local.h2))
##            rho     2.5 %    97.5 %
## [1,] 0.3245192 0.3007081 0.3479261
#dgn v. wright peripheral blood
ggplot(nona,aes(x=local.h2,y=h2PB))+geom_point(alpha=1/4)+geom_smooth()+theme_bw()

res<-cor.test(nona$local.h2,nona$h2PB,method='s')
CIrho(res$estimate,length(nona$local.h2))
##             rho      2.5 %     97.5 %
## [1,] 0.07345576 0.04716111 0.09964868
#wright v. price blood
ggplot(nona,aes(x=h2PB,y=h2bloodcis))+geom_point(alpha=1/4)+geom_smooth()+theme_bw()

res<-cor.test(nona$h2PB,nona$h2bloodcis,method='s')
CIrho(res$estimate,length(nona$local.h2))
##            rho      2.5 %    97.5 %
## [1,] 0.1064164 0.08025564 0.1324306

plot all 3 blood

nona <- all[complete.cases(all),]
dim(nona)
## [1] 5518   10
#dgn v. price blood, color by wright
ggplot(nona,aes(x=local.h2,y=h2bloodcis,col=h2PB))+geom_point()+geom_smooth()+theme_bw()+ scale_colour_gradient(low="pink",high="blue")

#dgn v. wright peripheral blood, color by price
ggplot(nona,aes(x=local.h2,y=h2PB,col=h2bloodcis))+geom_point()+geom_smooth()+theme_bw()+ scale_colour_gradient(low="pink",high="blue")

#wright v. price blood, color by dgn
ggplot(nona,aes(x=h2PB,y=h2bloodcis,col=local.h2))+geom_point()+geom_smooth()+theme_bw()+ scale_colour_gradient(low="pink",high="blue")

cor tables

round(cor(all[,-1],use='pairwise.complete.obs',method="s"),3)
##            local.h2 h2adip h2blood h2adipcis h2adiptra h2bloodcis
## local.h2      1.000  0.126   0.317     0.134    -0.006      0.307
## h2adip        0.126  1.000   0.335     0.328     0.617      0.231
## h2blood       0.317  0.335   1.000     0.246     0.140      0.392
## h2adipcis     0.134  0.328   0.246     1.000    -0.415      0.235
## h2adiptra    -0.006  0.617   0.140    -0.415     1.000      0.059
## h2bloodcis    0.307  0.231   0.392     0.235     0.059      1.000
## h2bloodtra    0.041  0.172   0.629     0.089     0.138     -0.339
## h2PB          0.083  0.111   0.190     0.066     0.074      0.124
## gtexWB        0.281  0.051   0.110     0.060    -0.007      0.120
##            h2bloodtra  h2PB gtexWB
## local.h2        0.041 0.083  0.281
## h2adip          0.172 0.111  0.051
## h2blood         0.629 0.190  0.110
## h2adipcis       0.089 0.066  0.060
## h2adiptra       0.138 0.074 -0.007
## h2bloodcis     -0.339 0.124  0.120
## h2bloodtra      1.000 0.102  0.004
## h2PB            0.102 1.000  0.042
## gtexWB          0.004 0.042  1.000
#Blood only
allblood <- dplyr::mutate(all,dgnWB=local.h2,Price=h2bloodcis,Wright=h2PB) %>% dplyr::select(dgnWB,gtexWB,Price,Wright)
#Pearson
round(cor(allblood,use='pairwise.complete.obs'),3)
##        dgnWB gtexWB Price Wright
## dgnWB  1.000  0.469 0.391  0.150
## gtexWB 0.469  1.000 0.249  0.118
## Price  0.391  0.249 1.000  0.156
## Wright 0.150  0.118 0.156  1.000
#Spearman
round(cor(allblood,use='pairwise.complete.obs',method="s"),3)
##        dgnWB gtexWB Price Wright
## dgnWB  1.000  0.281 0.307  0.083
## gtexWB 0.281  1.000 0.120  0.042
## Price  0.307  0.120 1.000  0.124
## Wright 0.083  0.042 0.124  1.000

cor tables - intersection of the 4 whole blood studies

dim(nona)
## [1] 5518   10
round(cor(nona[,-1],use='pairwise.complete.obs',method="s"),3)
##            local.h2 h2adip h2blood h2adipcis h2adiptra h2bloodcis
## local.h2      1.000  0.113   0.326     0.125    -0.012      0.325
## h2adip        0.113  1.000   0.272     0.374     0.572      0.234
## h2blood       0.326  0.272   1.000     0.261     0.065      0.450
## h2adipcis     0.125  0.374   0.261     1.000    -0.409      0.275
## h2adiptra    -0.012  0.572   0.065    -0.409     1.000      0.019
## h2bloodcis    0.325  0.234   0.450     0.275     0.019      1.000
## h2bloodtra    0.026  0.107   0.580     0.064     0.105     -0.337
## h2PB          0.073  0.050   0.147     0.041     0.041      0.106
## gtexWB        0.271  0.041   0.087     0.057    -0.023      0.130
##            h2bloodtra  h2PB gtexWB
## local.h2        0.026 0.073  0.271
## h2adip          0.107 0.050  0.041
## h2blood         0.580 0.147  0.087
## h2adipcis       0.064 0.041  0.057
## h2adiptra       0.105 0.041 -0.023
## h2bloodcis     -0.337 0.106  0.130
## h2bloodtra      1.000 0.066 -0.035
## h2PB            0.066 1.000  0.032
## gtexWB         -0.035 0.032  1.000
#Blood only
allblood <- dplyr::mutate(nona,dgnWB=local.h2,Price=h2bloodcis,Wright=h2PB) %>% dplyr::select(dgnWB,gtexWB,Price,Wright)
#Pearson
round(cor(allblood,use='pairwise.complete.obs'),3)
##        dgnWB gtexWB Price Wright
## dgnWB  1.000  0.461 0.407  0.134
## gtexWB 0.461  1.000 0.255  0.111
## Price  0.407  0.255 1.000  0.149
## Wright 0.134  0.111 0.149  1.000
#Spearman
round(cor(allblood,use='pairwise.complete.obs',method="s"),3)
##        dgnWB gtexWB Price Wright
## dgnWB  1.000  0.271 0.325  0.073
## gtexWB 0.271  1.000 0.130  0.032
## Price  0.325  0.130 1.000  0.106
## Wright 0.073  0.032 0.106  1.000